/*
 * InlineStream.cpp
 *
 *  Created on: Jan 5, 2011
 *      Author: blyth
 */

#include "InlineStream.h"
#include "PeakReference.h"
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/RandGauss.h"
#include <cmath>
#include <iostream>
using namespace std;

InlineStream::InlineStream(int length, int numPeaks, long seed, double* refPeak, int refSize) {
	this->length = length;
	this->numPeaks = numPeaks;
	randEngine = new CLHEP::MTwistEngine(seed);
	randGauss = new CLHEP::RandGauss(randEngine);
	this->refPeak = refPeak;
	this->refSize = refSize;
	stream = 0;
	peakInfo = 0;

	reRandomize();
}

InlineStream::~InlineStream() {
	if(stream)
		delete stream;
	if(peakInfo){
		for(int i = 0;i < numPeaks;i++)
			delete peakInfo[i];
		delete peakInfo;
	}
	delete randGauss;
}

void InlineStream::reRandomize(){
	if(stream)
		delete stream;
	if(peakInfo){
		for(int i = 0;i < numPeaks;i++)
			delete peakInfo[i];
		delete peakInfo;
	}

	stream = new int[length];
	for(int i = 0;i < length;i++)
		stream[i] = 0;

	peakInfo = new double*[numPeaks];
	for(int i = 0;i < numPeaks;i++){
		peakInfo[i] = new double[2];
		peakInfo[i][0] = randEngine->flat()*15;//*16384;
		peakInfo[i][1] = randEngine->flat()*(length-refSize);
	}

	PeakReference PR(refPeak, refSize);
	for(int i = 0;i < numPeaks;i++){
		int start;
		double offset = peakInfo[i][1]-(int)peakInfo[i][1];
		if(offset == 0)
			start = (int)peakInfo[i][1];
		else
			start = (int)peakInfo[i][1]+1;
		for(int j = 0;j < refSize && start+j < length;j++){
			stream[start+j] += PR.cubicInterpolationAt(j+offset)*peakInfo[i][0];
		}
	}

	addNoise();
}

void InlineStream::addNoise(){
	for(int i = 0;i < length;i++)
		stream[i] += (int)(randGauss->fire(0.,2.64575131));
}

void InlineStream::DoubleHitStream(int sizePeakOne,int sizePeakTwo,int seperation){
	if(stream)
		delete stream;
	if(peakInfo){
		delete peakInfo[1]; delete peakInfo[2]; delete peakInfo;
	}

	stream = new int[refSize*4+seperation];
	for(int i = 0 ; i<refSize*4+seperation ; i++)
		stream[i]=0;

	peakInfo = new double*[2];
	peakInfo[0] = new double[2];
	peakInfo[0][0] = sizePeakOne;
	peakInfo[0][1] = randEngine->flat()*(refSize);
	peakInfo[1] = new double[2];
	peakInfo[1][0] = sizePeakTwo;
	peakInfo[1][1] = peakInfo[0][1]+seperation;

	PeakReference PR(refPeak, refSize);
	for(int i = 0 ; i<2 ; i++){
		int start;
		double offset = peakInfo[i][1]-(int)peakInfo[i][1];
		if(offset == 0)
			start = (int)peakInfo[i][1];
		else
			start = (int)peakInfo[i][1]+1;
		for(int j = 0 ; j<refSize && start+j<refSize*4+seperation ; j++)
			stream[start+j] += PR.cubicInterpolationAt(j+offset)*peakInfo[i][0];
	}
}
